In MayaGans/CoreMicro: Testing Widely Used Core Microbiome Assignment Methods

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Load packages

library(CoreMicro)
library(patchwork)
library(phyloseq)
library(ggplot2)
library(tidyverse)

Create Human OTU and Arab OTU objects from data included in package

human_stool<-CoreMicro::human_stool
names(human_stool)[1] <- "X"
human_otu <- human_stool


#read in file
arabidopsis_rhizo<-CoreMicro::arabidopsis_rhizo
names(arabidopsis_rhizo)[1] <- "X"
arab_otu <- arabidopsis_rhizo

Figure 1

Define core_plots.core_methods function for plotting figure

core_plots.core_methods  <- function(combined_otu_data,
                                     low = "#EAE6f3",
                                     high = "#432976",
                                     legend_title = "Taxa Count") {

    combined_otu_data %>%
      ggplot() +
      aes(x = log(.data$Mean), y = .data$CV, color = as.factor(.data$value)) +
      geom_hex(bins = 30) +
      scale_color_manual(values = c("lightgray", "black")) +
      scale_fill_gradient(low = low, high = high) +
      theme_bw() +
      facet_grid(.~ .data$name, labeller = label_wrap_gen(multi_line = TRUE)) +
      guides(color = FALSE) +
      ylab("Coefficient of Variance") +
      labs(fill = legend_title)

}

Create panels of figure 1 for both datasets

dat <- CoreMicro::core_methods(human_otu)

levels(dat$name)[1] <- 'Abundance'
levels(dat$name)[2] <- 'Occupancy'
levels(dat$name)[3] <- 'Abundance and Occupancy'
levels(dat$name)[4] <- 'Hard Cutoffs'

p1.1 <-core_plots.core_methods(dat, legend_title = "HMP Taxa Count")  +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +                                                           
  theme(axis.text=element_text(size=20)) + 
  theme(axis.title=element_text(size=20)) + 
  theme(strip.text.x = element_text(size = 10))  + theme(legend.text=element_text(size=22)) + theme(legend.title = element_text(size = 24))




  dat <- core_methods(arab_otu)

levels(dat$name)[1] <- 'Abundance'
levels(dat$name)[2] <- 'Occupancy'
levels(dat$name)[3] <- 'Abundance and Occupancy'
levels(dat$name)[4] <- 'Hard Cutoffs'

 p1.2 <-  core_plots.core_methods(dat, high= "#1B301E", 
             low = "#E1E4E1", 
             legend_title = "Arabidopsis Taxa Count")+                                                           
  theme(axis.text=element_text(size=20))+ 
  theme(axis.title=element_text(size=20)) + 
  theme(strip.text.x = element_text(size = 10)) + theme(legend.text=element_text(size=22)) + theme(legend.title = element_text(size = 24))   

figure1 <- p1.1 / p1.2 #+ patchwork::plot_annotation(tag_levels = "A")

figure1